A  NON-HOMOGENEOUS  BINOMIAL  MODEL  FOR  THALAMIC 

OSCILLATIONS 


T.  Wang  and  J.  Muthuswamy 

Department  of  Bioengineering,  Arizona  State  University,  Tempe,  AZ  85287-9709 

Email:  jit@asu.edu 


ABSTRACT-  The  thalamic  ventral  posterior  lateral  neurons 
(VPL)  respond  to  somatosensory  stimulation  with  a  burst  of 
action  potentials  followed  by  a  periodic  oscillation  at  the 
spindle  frequency.  This  study  aims  to  build  a  statistical  model 
to  quantify  the  multi -unit  behavior  and  explain  putative 
underlying  mechanisms.  Multi -unit  data,  comprising  4  or  5 
different  neurons,  were  collected  from  anesthetized  adult  rats 
(n=2)  by  positioning  a  microelectrode  in  the  ventral  posterior 
lateral  (VPL)  nuclei  of  the  thalamus.  Using  an  observation 
window  of  1  ms  and  assuming  that  neuronal  firing  is 
uncorrelated  within  this  window,  the  firing  rate  of  the 
neurons  can  be  successfully  modeled  by  using  a  non- 
homogeneous  binomial  model  with  N=1  (with  99.5% 
confidence).  Using  maximum  likelihood  estimator  (MLE)  of 
the  parameter  p ,  statistically  consistent  prediction  of  the 
parameters  of  non-homogeneous  binomial  model  was  made 
using  a  minimum  of  50  stimulus -response  pairs.  The  inter¬ 
stimulus  interval  histograms  of  the  individual  neuronal  firing 
indicate  a  possible  stochastic  resonance  behavior  that  will 
model  the  spindles  in  thalamus.  Our  model  offers  a 
statistically  elegant  description  of  oscillations  in  neuronal 
action  potential  data  and  can  in  general,  be  used  to  track 
changes  in  the  neuronal  dynamics  with  function  or 
dysfunction. 
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I.  Introduction 

Neuronal  oscillations  have  been  observed  in  several 
regions  of  the  brain  either  spontaneously  or  in  response  to 
periodic  stimuli.  However,  quantification  of  these 
rhythmic  oscillations  has  been  confined  to  averaging 
histograms  using  auto  -correllograms  or  post-stimulus 
histograms  (for  stimulus-response  pairs)  [1].  These 
histograms  test  for  consistency  in  the  observation  but  give 
no  insights  into  the  generating  process,  nor  do  they  test  for 
statistical  significance.  An  earlier  effort  by  Konig  [2] 
attempts  to  fit  a  generalized  Gabor  function  to  histograms 
of  periodic  oscillations  using  least  squares.  The  parameters 
of  the  Gabor  function  then  quantify  the  oscillations. 
However,  this  model  still  does  not  describe  the  statistics  of 
the  generating  process  and  may  not  be  useful  for 
mechanistic  insights. 

This  study  focuses  on  multi-unit  oscillations  observed 
due  to  a  reciprocal  interaction  between  reticular  thalamic 
and  dorsal  thalamic  nuclei  in  response  to  somatosensory 
stimulation  [3].  Cortex-thalamus -cortex  loops  generate  a 
powerful  and  coherent  feedback  onto  the  thalamus, 
resulting  in  highly  coherent  oscillations  [4],  and  the 


synchrony  is  also  maintained  by  intrathalamic  connectivity 
[5].  We  investigate  appropriate  statistical  models  that  can 
describe  oscillations  recorded  from  dorsal  thalamic 
neurons  and  then  attempt  to  make  consistent  estimates  of 
the  descriptors. 

II.  Methodology 

Multi-unit  data  was  collected  from  adult  Wistar  rats 
(n=2).  Tungsten  electrodes  (FHC  Inc.,  Bowdoinham,  ME) 
were  implanted  in  the  thalamic  ventral  posterior  lateral 
(VPL)  area  for  recording  multi-unit  data.  The  rat  was 
under  barbiturate  anesthesia  (sodium  pentobarbital)  during 
the  experiment.  Periodic  somatosensory  stimulation  was 
given  every  3  sec.  Data  was  acquired  using  Spike2 
(Cambridge  Electronic  Devices,  Cambridge,  UK).  The 
units  were  sorted  using  a  template  matching  algorithm 
within  Spike2. 

Time -series  analyses  were  performed  on  the  successive 
spike -counts  measured  with  bin  width  of  1  msec.  The 
successive  firing  event  of  the  neuron  is  not  independent 
due  to  the  modulation  of  the  network  of  RT  (reticular- 
thalamic)  and  VPL  neurons,  so  the  Poisson  model  is  not 
suited  here.  But  within  a  small  time  window  of  1  msec, 
which  is  small  enough  for  just  one  action  potential,  we 
fitted  a  non-homogeneous  binomial  model.  A  Maximum 
Likelihood  Estimator  (MLE)  is  used  here  to  estimate  the 
probability  of  the  firing  P(t),  which  changes  with  time 
after  the  stimulus.  Wilcoxon  rank  sum  test  is  used  to 
evaluate  the  statistical  hypothesis.  A  p-\ alue  (0<£><1) 
much  larger  than  zero  indicates  a  good  match  between  the 
model  and  the  data.  The  evolution  of  estimated  mean  and 
variance  of  the  firing  probability  was  investigated  every 
five  stimulus  to  determine  the  minimum  number  of 
stimulus-response  pairs  needed  to  make  a  statistically 
consistent  estimate  of  parameter  P(t). 

Cross -correlation  and  power  spectrum  is  used  to 
investigate  the  relationship  of  P(t)  of  two  different 
neurons.  To  investigate  the  burst  behavior  in  each  spindle, 
we  separate  the  3  sec  inter-stimulus  interval  into  30 
segments  and  plot  the  ISI  of  each  segment.  The  ISI  plots 
were  also  examined  for  the  presence  of  stochastic 
resonance. 
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III.  Results 


Thalamic  oscillations  have  two  components,  a 
modulating  spindle  oscillation  which  is  caused  by  cortex- 


Fig.  1.  Post-Stimulus  Histogram  (PSH)  of  (a). 
Thalamic  oscillation  (b).  inner  spindle  burst 


thalamus -cortex  loops  with  frequency  of  about  10  Hz  (Fig. 
1(a)),  and  the  inner  spindle  bursts  (Fig.  1(b). 

Using  one  milli-second  time  bin  with  at  most  one  spike 
(action  potential)  in  each  window,  the  whole  oscillation 
can  be  modeled  as  a  non-homogenous  binomial  process, 
where  P(t)  (the  probability  of  firing  at  time  t)  changes  with 
time  (Fig.  1(a)).  The  ^-values  from  the  Wilcoxon 
hypothesis  tests  at  different  times  after  stimulation  are 
shown  in  Fig.  2.  The  large  p-values  during  the  times  of 
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Fig.  2.  Wilcoxon  hypothesis  test  result 

action  potential  firing  indicate  a  good  fit  to  the  non- 
homogeneous  binomial  model. 

Using  a  maximum  likelihood  estimator  (MLE),  we 
estimated  the  evolution  of  mean  (P(t))  and  variance  of  the 
firing  rate  at  specific  time -points  after  stimulus.  The 
estimates  are  plotted  against  increasing  numbers  of 
stimulus-response  pairs  in  Figs.  3(a)  and  (b).  The  two  plots 
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Fig.  3.  Evolution  of  meari^nd  variance  at  (a)  the 
first  peak  (n=16  ms)  and  (b)  the  second 
peak(n=122  ms) 

shown  correspond  to  16  msec  and  122  msec  after  stimulus 
(corresponding  to  the  first  and  second  peak  in  Fig.  1(a)). 
Fig.  3  shows  that  the  neuronal  firing  is  a  stationary  random 
process.  With  approximately  50  stimulus-response  pairs 
we  can  make  a  good  prediction  of  the  mean  ( P(t )). 

After  sorting  the  multi-unit  data,  we  investigate  the 
dynamics  of  each  neuron.  The  power  spectra  of  P(t)  of  two 
neurons  are  shown  in  Fig.  4  with  very  similar  profiles  and 
fundamental  frequencies.  The  cross  correlation  of  P(t) 
between  the  two  neurons  indicates  a  1ms  phase  shift.  Figs. 
5(a)  and  5(b)  indicate  that  the  most  frequent  occurrence  of 
the  inter-spike  interval  is  around  3  msec  for  both  neurons. 
There  are  also  non-periodic  fluctuations  at  4ms,  5ms,  and 
6ms. 


Fig.  4.  Power  spectrum  of  P(t)  of  two  neurons 


IV.DlSCUSSION  AND  CONCLUSIONS 


A  non-homogeneous  binomial  model  has  been  used  to 
describe  the  statistical  behavior  of  the  thalamic 
oscillations.  Consistent  estimates  of  the  probability  of 
firing,  P(t)  can  be  obtained  from  a  minimum  of  50 
stimulus-response  pairs  by  using  a  maximum  likelihood 
estimator  (MLE).  This  provides  an  elegant  statistical 
description  of  the  thalamic  response  as  whole  and  neuronal 
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Figure  5:  ISI  distribution  of  each  spindle  after  stimulation 
(a),  neuron  1.  (b)  neuron  2.  The  spindle  numbers 
correspond  to  the  peaks  in  the  histogram  following  the 
stimulus.  Each  peak  corresponds  to  a  spindle.  The  time 
labels  corresponds  to  the  time  after  stimulation. 

oscillations  in  particular.  Conventional  methods  to 
estimate  periodicities  in  neuronal  oscillations  involve 
averaging  of  approximately  100  response  histograms.  The 
only  known  effort  to  quantify  neuronal  oscillations 
involves  fitting  a  known  periodic  function  to  the  averaged 
histogram  [2],  which  does  not  describe  the  statistics  of  the 
generating  process  . 

Results  from  cross -correlation  indicate  that  the  neurons 
fire  synchronously  with  a  constant  time  delay.  In  the  ISI, 
there  is  a  strong  first  peak  at  the  time  of  about  3msec, 
which  means  the  neuron  most  probably  wait  about  3msec 


to  fire  again.  It  implies  that  in  VPL  there  should  exist  an 
intrinsic  periodic  driving  force,  but  it  is  not  strong  enough 
and  must  have  the  help  of  intrinsic  noise  to  make  the 
neuron  fire  [6].  In  fact,  Na+  channels  provide  a  possible 
source  of  intrinsic  noise  [7],  and  the  network  of  oscillators 
may  provide  the  intrinsic  periodic  driving  force  for  each 
neuron  [8].  Future  work  needs  to  be  done  by  using 
stochastic  resonance  to  understand  the  underlining 
mechanism  of  the  non-homogenous  binomial  model. 
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